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We propose a new formulation of the fluctuating lattice Boltzmann equation that is consistent 
with both equilibrium statististical mechanics and fluctuating hydrodynamics. The formalism is 
based on a generalized lattice-gas model, with each velocity direction occupied by many particles. 
We show that the most probable state of this model corresponds to the usual equilibrium distribution 
of the lattice Boltzmann equation. Thermal fluctuations about this equilibrium are controlled by 
the mean number of particles at a lattice site. Stochastic collision rules are described by a Monte 
Carlo process satisfying detailed balance. This allows for a straightforward derivation of discrete 
Langevin equations for the fluctuating modes. It is shown that all non-conserved modes should be 
thermalized, as first pointed out by Adhikari et al; any other choice violates the condition of detailed 
balance. A Chapman-Enskog analysis is used to derive the equations of fluctuating hydrodynamics 
on large length and time scales; the level of fluctuations is shown to be thermodynamically consistent 
with the equation of state of an isothermal, ideal gas. We believe this formalism will be useful in 
developing new algorithms for thermal and multiphase flows. 



I. INTRODUCTION 

Lattice Boltzmann (LB) methods [l], [2[ have become a 
popular tool for simulating hydrodynamics, particularly 
in complex geometries. The underlying model is a regular 
lattice of sites r, combined with a small set of velocity 
vectors which, within one time step h, connect a given 
site with some of its neighbors. The set of velocities 
is chosen to be compatible with the symmetry of the 
lattice. The basic dynamical variables are real-valued 
populations ti^; in the present paper, we will consider 
as the mass density associated with the velocity c^. The 
LB algorithm is then described by the update rule 

n,{f+c,h,t + h)^n*{r,t)^n,{f,t) + A,{n,{r,t)}, (1) 

where {rii} denotes the complete set of populations. The 
{ni{f, t)} at each site are first re-arranged in a "collision" 
step, described by A^, and then propagated along their 
respective links. The hydrodynamic fields, mass density 



i 

and momentum density 



(2) 



(3) 



are moments of the discrete velocity distribution Ui (r, t) , 
while the fluid velocity is given by 

u{r,t)=j{r,t)/p{r,t). (4) 

The collisions conserve mass and momentum, hence 

^ A, = ^ A,c, - 0. (5) 



The algorithm thus satisfies important requirements for 
simulating hydrodynamic flows — mass and momentum 
conservation, and locality — but lacks Galilean invari- 
ance due to the finite number of velocities. Full rotational 
symmetry is also lost, but by a suitable choice of veloc- 
ity set, isotropic momentum transport can be recovered 
on sufficiently large (hydrodynamic) length scales. Nev- 
ertheless, the finite number of velocities always confines 
the method to flows with small Mach number u/cg ^ 1. 
The speed of sound Cs is of order b/h, where b is the 
lattice spacing, or of order \ci\. 

Most of the LB literature deals with deterministic col- 
lision rules, with Aj describing a linear relaxation of the 
distribution {rii} towards the local equilibrium [3, 



(p, u) = pa"^ I + 
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u 
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(6) 



where a'^' > is the weight associated with the speed \ci\. 
The viscosity of the LB fluid is controlled by the choice 
of relaxation rates. 

However, to simulate Brownian motion of suspended 
particles, thermal fluctuations must be included. At 
the hydrodynamic level, this means adding uncorrelated 
noise to the fluid stress tensor In Refs. [1, 0, @] an 
analagous fluctuating LB model was introduced by mak- 
ing Ai a stochastic variable, but in such a way that the 
noise was only applied to the modes (linear combinations 
of {rii}) related to the viscous stress tensor 



nneq \ ^ nen 



here a and /3 denote Cartesian components and n 



(7) 
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is the non-equilibrium distribution. Although 



this procedure is correct in the hydrodynamic limit 0, [ 
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it provides poor thermalization on smaller length scales, 
as was first observed by Adhikari et al. [10t\ . They in- 
troduced a thermalization procedure which applies to 
all non-conserved modes, with significantly improved nu- 
merical behavior at short scales [fO]. The procedure was 
derived by considering a fluctuating LB model, making 
explicit use of the transformation between the popula- 
tions {n.;} and the modes 

The purpose of the present paper is to re-derive the 
stochastic updating rule of Ref. [f 0] from a generalized 
lattice-gas model. The novelty of our formulation lies in 
the introduction of an ensemble of population densities 
at each grid point, so that a fluctuating LB simulation is 
a single realization of this ensemble. There follows nat- 
urally a probability distribution, P({nj}), for the set of 
populations {n.i} at a position r and time t. The equilib- 
rium distribution at a single site can be derived by max- 
imizing P subject to the constraints of fixed mass and 
momentum densities, p and j. This distribution agrees 
with the standard equilibrium distribution for LB models 
[Eq. ^] up to terms of order u^. A similar procedure 
has been followed in deriving i7-theorems for LB mod- 
els [13, [13] J but these papers were not concerned with 
fluctuations. 

A coarse-graining of the microscopic collision opera- 
tor leads to a Langevin description for the non-conserved 
degrees of freedom. However these stochastic collisions 
may also be viewed as a Monte Carlo procedure p^ . 
satisfying the principle of detailed balance governed by 
P {{ui}). The procedure of Refs. 0, Q can be shown to 
violate detailed balance, while the improved version of 
Ref. satisfies it. 

In summary, our goal is to reconnect the lattice Boltz- 
mann equation with its lattice gas origins, and thus 
to establish a firm statistical mechanical foundation for 
stochastic LB simulations, as well as the usual connec- 
tion to fluctuating hydrodynamics 0, • We believe this 
provides a comparable theoretical framework to that al- 
ready available for other stochastic simulation methods, 
such as dissipative particle dynamics [l^ and stochas- 
tic rotation dynamics This formulation also offers 
the possibility for future modifications and generaliza- 
tions, for example to thermal flows [l8|, or models with 
nonideal equations of state (lol . [20| , or multi-component 
mixtures [2]| . 

The paper is organized as follows: In Sec. |TT] we de- 
scribe the underlying lattice-gas model, derive the proba- 
bility distribution P ({rii}), and show that the most prob- 
able value for {n;} is equivalent to Eq. In Sec. lIIII we 
consider small fluctuations around the equilibrium distri- 
bution. We show they are approximately Gaussian dis- 
tributed, with the level of thermal fluctuations governed 
by the degree of coarse-graining: a given amount of mass 
on a lattice site can be distributed between many parti- 
cles, in which case the fluctuations are small, or between 
few, in which case they are large. In this way we can ad- 
just the level of fluctuations, while keeping the tempera- 
ture fixed. In Sec. IIVI we construct a stochastic collision 



operator such that detailed balance is satisfied. From 
this, we derive the stochastic stresses at an individual 
site. In Sec. |V] we apply the Chapman-Enskog proce- 
dure to the algorithm in order to find the behavior on 
the hydrodynamic scale; the deterministic and stochas- 
tic terms are here treated on an equal basis 22]. We 
then find that, on the macroscopic scale, the procedure 
yields exactly the stress correlations given by Landau and 
Lifshitz [3] . Section |Vl] then discusses how to choose pa- 
rameters for a coupled particle-fluid system. Section IVIII 
summarizes our conclusions. 



II. SINGLE-SITE PROBABILITY 
DISTRIBUTION 

Historically, the lattice-Boltzmann model [1, [2^ de- 
veloped from earlier work on lattice-gas (LG) mod- 
els [23, UHl , in which each velocity direction was occupied 
by at most one particle. We imagine a generalized lattice- 
gas model (GLG) where each velocity direction can be 
occupied by many particles. Each particle has the same 
mass, but different velocity directions may have different 
mean populations, even in a fluid at rest. The micro- 
scopic state of the system at any given site is specified 
by a set of integers {i^i} giving the occupancies of each 
direction. Then the update of the GLG is analagous to 
the standard LG or LB models, but with an integer as 
opposed to a Boolean or real variable: 

iy,{r + Cih,t + h) = u*{r,t) = v,{r,t) + Ki{ui{r,t)), (8) 

where A.^ operates on {vi\ to compute the change in pop- 
ulation v* — Vi. While collisions may be both determin- 
istic and microscopically reversible, we shall assume only 
that the collision operator satisfies detailed balance. 

Without considering the collision rules in detail, we 
construct an equilibrium distribution from the the fol- 
lowing thought experiment. Consider a "velocity bin" 
I, related to one particular site r. This bin is placed in 
contact with a large reservoir of particles, such that the 
number of particles in the bin, vi, is a random variable. 
The probability for a particle to be in the reservoir is 
close to unity, and the probablility to be in the bin is 
small. Therefore, Vi follows a Poisson distribution, with 
a mean number of particles vi, 

P(^,) = ^e-'^-, (9) 

and a variance 

{u^) {u,f = u,. (10) 

Let TOp be the mass of a particle and /i = rup /¥, with d 
the spatial dimension of the system. Then Ui — fivi, and 
hence 

{n^)-{n,f ^^^(n,). (11) 
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The fluctuations in mass density at a site are controlled 
by the mass of an LB particle: small nip means that the 
mass is distributed onto many particles, and therefore 
fluctuations are small. For fixed rup, /i (and therefore 
the level of fluctuations) becomes large as b decreases. 
This is natural, since a fine spatial resolution means fewer 
particles per cell, and larger fiuctuations relative to the 
mean. 

If we now imagine sampling each velocity with an in- 
dependent reservoir, but taking only those sets of popu- 
lations which produce specific values for the total mass 
and momentum, the probability density for the occupa- 
tion numbers is (except for normalization) 



(12) 



Using Stirling's approximation for 3> 1, we can write 
the distribution in terms of the entropy associated with 
the occupation numbers, 

S [{vi}) = - ^ {vi In Vi - Vi - Vi In Vi + Vi) , (13) 



and the constraints: 

P{{vi}) cx e,qy[S {{vi})] 



(14) 



The equilibrium distribution, v'l'' , can be found by maxi- 
mizing S, treating Vi as a continuous variable, and taking 
into account the mass and momentum constraints via La- 
grange multipliers, Ap and Aj respectively: 



dS_ 

dvi 



Ap + Aj • Ci = 0, 



0, 


(15) 


0, 


(16) 


0. 


(17) 



It should be noted that this procedure is closely related 
to the determination of an entropy function for the LB 
equation [l3|. Equation can be solved to give the 
equilibrium populations in terms of the Lagrange multi- 
pliers. 



vl"^ = Vi exp {\p + A J • c,;^ , 



(18) 



which are then determined from the constraints, 
Eqs. (fTB)) and p7|) . subsituting v^'^ for Vi. 

The mean populations in the absence of constraints, 
{vi}, can be expressed in terms of the mean number of 
particles at a site. 



Vi = va ^ , 



(19) 



where v = Vi. The symmetry of the lattice constrains 
the weights, a*^', to be dependent on the speed of the 
particle, but not its direction. Thus for a lattice with 
cubic symmetry. 
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i 
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(20) 
(21) 
(22) 
(23) 



where Sap is the Kronecker delta, and 172 is a constant 
with units (b/h)'^. 

We seek an approximate expression for the equilibrium 
distribution in the limit that Aj-Ci is small 14| . To second 
order in Aj, the mass and momentum constraints yield: 



//Pe-^" 1 



(24) 
(25) 



Inserting these results into Eq. p8)) . we find the equilib- 
rium distribution can be written in the form of Eq. 



^ve "02^^^ pu. 



pa^^ 1 




(26) 



For the sake of completeness, we now briefly mention 
the well-known procedure [1, Q to determine the weights 
a"^' such that the LB model is consistent with hydrody- 
namics. This requires that the second moment of the 
equilibrium distribution. 



na/3 ^^nl^^C^aCip 



(27) 



should equal the Euler stress p5ap + pUaUp, with the 
pressure given by the ideal gas equation of state, p = 
pksT / Trip, where ks is Boltzmann's constant and T is the 
absolute temperature. For an isothermal gas of particles 
of mass nip, ksT = nipC^, and therefore the equation of 
state is also given hy p = pc^, with Cs the speed of sound. 

To evaluate 11^'^ we require the fourth moment of 
{a'^'}, which from cubic symmetry must be of the form 

y^a'^'CiaCi^Ci^CiS = KiSafj-fS (28) 
i 

where Safjjs is unity if all four indexes are the same and 
zero otherwise; K4 and (74 have units of (b/h)'^. Con- 
sistency between Eq. (I?7l) and the Euler stress requires 
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that: 



0-2 

(74 
K4 



'2> 



= 0. 



(29) 
(30) 
(31) 



These conditions, together with the normaHzation con- 
dition, ~ 1' determine the weights uniquely for a 
model with three different speeds. For example, for the 
D3Q19 model [| (19 velocit ies on a three-dimensional 



simple cubic lattice). 



1/3 for the stationary parti- 



cles, a = 1/18 for the six nearest -neighbor directions, 
and = 1/36 for the twelve next-nearest neighbor 
directions: the sound speed is then cl = (l/3)(&//i)^. 



vectors [lO|, |26| , for which the kinetic or "ghost" modes 
have no projection on the equilibrium distribution. We 
consider only the non- equilibrium distribution, which we 
can write as an orthonormal transformation of the scaled 
variables, Xi: 



ruk 



E 



k 



(36) 
(37) 



where rrij is the amplitude of the jth mode, and the basis 
vectors satisfy the orthonormality conditions 



E 



^ki^li 



^ki- 



rn 



III. SINGLE-SITE FLUCTUATIONS 

We now consider the distribution of small fluctuations 
in the mass densities associated with each velocity direc- 
tion, nf'^'' = rii — nl'^ . Using the results of Appendix YK\ 
to incorporate the constraints, and converting from fluc- 
tuations in Vi to fluctuations in m, 



P{{nr}) (X exp 



2/in 



eq 



(32) 



The variance in the fluctuations depends on direction, 
but, since n^'^'^ is already a small quantity in comparison 
with rtl*, we will approximate the variance by the low- 
velocity limit. 



n-^ ^ pa '. 



(33) 



We now introduce normalized fluctuations Xi , via the def- 
inition 



n^'"^ — y/Jipa^Xi, (34) 
and transform Eq. ((32)l to the simplified expression 



P({x,}) oc exp l^-i^x,'^ 



(35) 



HZ 



It should be noted that the basis vectors e^i are differ- 
ent from the Cki defined in Ref. [1^ , since there the trans- 
formation was for unsealed variables, rij, rather than the 
scaled variables, Xi, used here. The essential underlying 
physics of the transformation is however unchanged; the 
present expressions are just a re-parametrization. The 
basis vectors are related to the weighted basis vectors 
used in Ref. 1261: 



^ki 



where Wk is the length of the k'th basis vector. 



Wk 



E' 



(39) 



(40) 



The hydrodynamic modes, mass density, momen- 
tum density, and stress, can be written in a model- 
independent form. Explicitly, 



for the mass mode, and 



1, 



(41) 



(42) 



for the momentum modes. Note that in our formalism 
Too and TOq {a — 1, . . . , d) are zero. 

In addition to the conserved modes, there are d{d + 
l)/2 viscous modes: one bulk mode, d—1 shear modes 
involving diagonal elements of the CiCi tensor, and d{d — 
l)/2 off-diagonal elements. The bulk stress mode is given 

by 



Eqs. ([6]), ((34|) . and ((35|) define the statistics of our fluc- 
tuating LB model. 

The LB collision operator can be conveniently repre- 
sented in terms of modes, which are linear combinations 
of the mass densities, {rii} [ll|. The basis vectors are 
constructed from orthogonal polynomials in the veloc- 
ity set {ci}. There is more than one possible choice for 
these basis vectors [2^ , and we use the "weighted" basis 



(43) 



where orthogonality to the mass mode is assured by 
Schmidt orthogonalization. There is a shear mode of the 
form 



1 



c2 V 2d{d - 1) 



(44) 
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TABLE L Basis vectors of the D2Q9 model. Each row cor- 
responds to a different basis vector, with the actual polyno- 
mial in Cia shown in the second column; the components of 
Cia = Ciah/b are normalized to unity. The orthonormal basis 
vectors iki can be obtained from the table using Eq. (|39|l . 
Cfei = •\/a'=i /wkeki- the normalizing factor for each basis vec- 
tor is in the third column. 



k 




■Wk 





1 


1 


1 




1/3 


2 


Ciy 


1/3 


3 


3cf -2 


4 


4 




4/9 


5 




1/9 


6 


(3cf - 4)c^x 


2/3 


7 


(3c? - 4)ci,, 


2/3 


8 


9c- - 15c? + 2 
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and d — 2 shear modes of the form {d > 2) 



2c2 



(' 



(45) 



together with additional modes formed by cychc permu- 
tations of the Cartesian indexes. The d{d — l)/2 off- 
diagonal shear stresses are of the form 



2 ^ix^iy 



(46) 



together with cyclic permutations. 

All these vectors are mutually orthogonal. Further or- 
thogonal vectors, whose span are the so-called kinetic or 
"ghost" modes, may be constructed in terms of higher- 
order polynomials of Ci [ll|; these are model dependent. 
A complete set of basis vectors for the D2Q9 and D3Q19 
LB models ^ are given in Tables U and |TT] respectively. 

Equation (|35p can be rewritten using Eqs. (|37p and 
([55)1 , to give the non-equilibrium probability distribution 
of the modes rrik, 

P({mfc}) cx exp I -i^m^. I ]^5(m,) 

\ k / i<d 

k>d / 



CX exp 



(47) 



There is no contribution to P from the conserved modes. 



IV. STOCHASTIC COLLISIONS AS A MONTE 
CARLO PROCESS 

In this section we construct a stochastic collision oper- 
ator, viewed as a Monte Carlo process, and consider the 



TABLE II; Basis vectors of the D3Q19 model. Each row 
corresponds to a different basis vector, with the actual poly- 
nomial in Cia shown in the second column; the components of 
Cia ~ Ciah/b are normalized to unity. The orthonormal basis 
vectors Cki can be obtained from the table using Eq. (I39|l . 
Cki = ^/a,'^' /wkGki- the normalizing factor for each basis vec- 
tor is in the third column. 
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1 




1/3 
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Ciy 


1/3 


3 




1/3 
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c?-l 


2/3 


5 


3c^ -c^ 


4/3 


6 


^ly ^iz 


4/9 


7 


C-ix Ciy 


1/9 


8 


Ciy Ciz 


1/9 


9 


Ciz Cix 


1/9 


10 


{3Ci - 5)cix 


2/3 


11 




2/3 


12 


(3Ci - 5)ciz 


2/3 


13 


{C-iy ~ C^z)Cix 


2/9 


14 


{Ciz Cj^x)Ciy 


2/9 


15 


{.Cix Ciy)Ciz 


2/9 


16 


3c| - 6c? -1- 1 


2 


17 


(2c?-3)(3c?,-c?) 


4/3 


18 


(2Cj 3)(C|y — Ciz) 


4/9 



local dynamics at the level of a single lattice site. In the 
next section (Sec. IV]) we will consider the global dynam- 
ics, through a Chapman-Enskog expansion. A determin- 
istic collision operator at the microscopic level is quite 
complicated to constructjCven for the simplest three- 
dimensional LG models j23|, and cannot be easily ex- 
tended to the larger number of particles in Eq. ([5]). Col- 
lision rules are much easier to construct at the Boltzmann 
level 0] ; the stochastic update from pre-coUision to post- 
collision populations, ni n*, is facilitated by making 
the transition between modes, m\, since each de- 

gree of freedom is then independent. Denoting a tran- 
sition probability between the pre- and post-collision 
states of a particular mode m by uj[m — > m*), the condi- 
tion of detailed balance, governed by the distribution in 
Eq. (gT]), reads 



ujirn —> m*) exp (-m*^/2) 
uj{m* m) exp {—m?/2) 



(48) 



A simulation at the hydrodynamic level does not need 
to satisfy this condition, and typically does not, but it 
is essential for a proper thermal equilibrium of the LB 
fluid. 

There are many possible realizations of Eq. : one 
well-known example is the Metropolis method, involving 
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a trial move followed by a stochastic acceptance or rejec- 
tion step to enforce detailed balance. Here we consider 
the linear relaxation model typically used in LB simula- 
tions, balanced by Gaussian noise: 

m* = 7m -|- Lpr, (49) 

where 7 is related to an eigenvalue of the linearized col- 
lision operator, 7 = 1 -I- A (see Eq. 8 of Ref. [2^), and r 
is a Gaussian random number with zero mean and unit 
variance. The dissipation parameter 7 is restricted by 
the linear stability limit, I7I < 1, with the case 7 < 
corresponding to "over-relaxation". Equation ((49)) has 
the technical advantage of being rejection-free, and the 
conceptual advantage of enabling an analytic calculation 
to be made at the Ghapman-Enskog level (see Sec. W\ . 

The parameter </? must be adjusted to satisfy de- 
tailed balance, Eq. (jiS]) . using the relation [Eq. (jiO]) ] 
r — ip^^{m* — 7m). Since the transition probability for 
TO TO* is identical to the probability for generating the 
value of r that gives to* from m, 

w(to ^ m*) = (27r(/?^) exp [— (m* — 7to)^/2(p^] . 

.(50) 

There is a similar expression for the reverse transition, 
TO* m, with TO* and to interchanged. From Eq. (|48ll . 
we then find that detailed balance is satisfied for 

<^=(l-7^)i/^ (51) 

Thus the case 7=1 corresponds to a conserved mode, 
while 7 = corresponds to m* being entirely random, 
with no memory of its previous value. 

Each mode, toj,, in the LB model is assigned its own 
relaxation rate 7^, subject to the constraints of symme- 
try and conservation laws; the conserved modes {k < d) 
require that 7^ = 1. For the bulk stress we choose a 
value 76, and for the {d + 2)(rf — l)/2 shear stresses a 
single value 7s. In Refs. 0, H, [iP'l the kinetic modes were 
updated with 7^ = 0, but it is possible to achieve more 
accurate boundary conditions with a propert tuning of 
the kinetic eigenvalues [1^, [2^ . Equation ([5T|) ensures 
that detailed balance is satisfied for all choices of jk- 
A purely deterministic LB model is obtained by setting 
ipk = for all modes; physically, this corresponds to the 
limit of nip ^ 0, or ^ 00. 

The original formulation of the fluctuating LB 
model 0, Q is obtained by setting jk = fk = for 
all the kinetic modes, but choosing the variance of the 
stresses according to Eq. ([?T|) . The kinetic modes are 
projected out at every time step by this collision rule, and 
Lu (mk — > TO^. = 0) = 1. However, there is no route back 
to the pre-coUisional state, lo (to^ = ^ to^) — 0, and 
detailed balance [Eq. is clearly violated. Neverthe- 
less, this model still yields the correct fluctuating hydro- 
dynamics in the limit of large length scales Q , as is shown 
by the analysis in Sec.|Vl Treating all the non-conserved 
modes on an equal basis [Tj satisfies detailed balance on 
all scales, and is entirely equivalent to Eqs. (|151) - ([5T|) . 



As a general rule, proper thermalization requires as 
many random variables as there are degrees of freedom 
in the system (not counting the conserved variables). 
When 7A: = for all kinetic modes, the deterministic LB 
model can be propagated forward in time knowing only 
the mass, momentum, and stress tensor at each lattice 
site. Thus it might appear that the {n^} play the role 
of auxiliary variables, and in this case the LB model has 
fewer degrees of freedom than when 7^ 7^ 0. However 
this is incorrect, since the deterministic dynamics with 
any 7^ = is pathological, in the sense that it is now im- 
possible to reconstruct the trajectory backwards in time 
(in contrast to the case when all 7^ ^0). Thus the num- 
ber of degrees of freedom is well-defined only in the case 
where all 7fe ^ 0. The ill-defined case (some 7^ = 0) is 
however a limiting case of the well-defined one, and thus 
continuity tells us that the number of random variables 
should be the same in both cases. 

The update rule in Eq. (|49p. with 7 > Oj_is an exact 
solution of a continuous Langevin equation |29l . [30| , 

^m = -Tm + £, (52) 
at 

with {i{t)) = and {i{t)i{t')) = 2T5{t - t'). Integrating 
Eq. (l52|) from t = to t = h [i. e. one LB time step) 
gives Eq. with 7 — exp(— F/i). The standard first- 

order Euler approximation to Eq. (j52p corresponds to 
7 1 — r/i, and is only valid for small Th. By contrast 
Eq. does not impose any restriction on the time step. 

For the Chapman-Enskog analysis in Sec. |Vl we will 
need the coUisional update of the non-equilibrium stress 
tensor, H^^'' = X^i equilibrium part of 

the stress is unchanged by the collision process. We first 
decompose H^^'' into a multiple of the unit tensor (bulk 
stress), and a traceless part (shear stresses), denoted by 
an overbar: 

nr/ = nr/ + ^n-'^„/3, (53) 

where we have used the Einstein summation convention 
for the Gartesian components. The change in the non- 
equilibrium stress tensor at a lattice site, due to collisions, 
can be determined from Eqs. (^5]) and ([5T|) . 

KT = isK7 + Ro.0. (54) 
n*r = 76nr + (55) 

The variables Raf3 are Gaussian random variables with 
zero mean; in addition R^p is traceless. The covariance 
matrix {RapR-^s) is determined by the variances of the 
stochastic stress modes. The calculation can be simpli- 
fied by observing that the matrix is a fourth rank tensor 
and is therefore isotropic by the symmetries of the LB 
model, 

{RapRfS) = Rl {SajSpS + SasSpj) ■ 

+ R2SapSyS, (56) 
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The unknown constants, Ri and i?2, can be determined 
from special cases. For example, in the D3Q19 model 
defined in Table El 

nsr = VMPc>7, (57) 

and therefore, from Eq. 

(Rly) ^ fipct {1 ^ j',) = Ri, (58) 
where the final equality follows from Eq. (|56p . Similarly, 

n^r-n"r = 2VMPc'"^6, (59) 

and therefore 

{{Ryy - Rzzf) - ^f^pct (1 - 7') = 4i?i, (60) 



This result is consistent with Eq. (|56p. which demon- 
strates that the fluctuating stresses are indeed isotropic. 
Finally, the fluctuations in the trace, Raa, are related to 
lb- 



Multiplying this equation by one of the basis vectors and 
summing over all the directions, we obtain the equations 
for the dynamics of the fluctuating LB model on the ti 
time scale, 

dti^n'l'^eki + dl^nl'^Caeki = h^^ ^ A^Cfei. (65) 



Note that we use the Cki basis vectors here, in conjunction 
with and {n^'^'^}, not the normalized basis vectors 

Cfci, which are for the {xi\. 

When applied to the conserved degrees of freedom, k < 
d, Eq. leads to the inviscid fluid equations: 

dt^P + dijo. = 0, (66) 
dtija + djs {pclSap + pUaUp) = 0, (67) 

Similarly, for the stress modes, d < k < (dP + 3(i)/2, we 
find: 



(68) 



The general expression for the covariance in the ran- 
dom stresses is 



{Ra/nRfS) 



= (1 - 7^ ) {Sa-ySps + SasSp-y) (62) 

2 

+ ^ (7s - lb) Sa^S^S- 



Evaluating the time derivatives in Eq. (|68p gives a simpli- 
fied expression for the non-equilibrium stress, apart from 
small terms of order Q, 



(69) 



This covariance matrix is different from the global fluc- 
tuations in stress, which are superposed onto the hydro- 
dynamic modes (Sec. \V\. 



V. CHAPMAN-ENSKOG EXPANSION 

In order to determine the behavior on hydrodynamic 
length and time scales, we apply the Chapman-Enskog 
method to the stochastic dynamics of the fluctuating LB 
model. We modify the derivation of Ref. @ to include 
thermal fluctuations: for an alternative procedure, see 
Ref. [Slj. Here, the expansion parameter e is used to sep- 
arate the lattice scale, r, from the hydrodynamic scale, 
ri = er. Thus da = e9^, with the notation da = d/dva- 

Since the collision operator is local in space and time, 
the non-equilibrium distribution is also taken to be of 
order e: n"^'^ — en], with nj of order unity in the e 
expansion, 



Thus, on the ti time scale, the viscous stresses fluctu- 
ate around a mean value that is slaved to the velocity 
gradient, 9^M/3 -I- d^Ua- 

The kinetic modes fluctuate around zero on the ti time 
scale, with at most a small correction of order u^: 



ml = rrik + O(m^). 



(70) 



The equilibrium distribution contains polynomials in Ci 
up to second order, and is thus automatically orthog- 
onal to the kinetic modes, which are made up of 3rd- 
order and 4th-order polynomials in ci . Since the equilib- 
rium distribution has no projection on the kinetic modes, 
the time-derivative in Eq. (|65p vanishes identically for 
k > (d^ -I- 3d)/2. However, the gradient term in Eq. ([65]) 
includes an additional factor of Ci : thus third-order poly- 
nomials survive, making small equilibrium contributions 
of order to the dynamics. 

At the order e^, the Boltzmann equation is 



■ en. 



(63) 



neq 



1 neq 



We use the usual multiple time scale expansion [32l |. 
dt = edti + e^dt2, to separate the convective (ti) and 
diffusive (^2) relaxation processes. The left-hand side of 
Eq. IT]) is expanded about (r , t) in a Taylor series with 
respect to h: to first order in e. 



Ciadi) n\ 



(64) 



+ 2^*1 (^*i + 

+ ^dl{dt,+c,adi)n'r''c,p = 0, (71) 

where the terms have been grouped to suggest the most 
expedient means of calculation. Since only the hydro- 
dynamic modes survive to the t2 timescale, we consider 
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just the modes up to fc = d. It follows immediately from 
Eq. ([7T|) and the conservations laws [Eas. ([S51) and ([57)1 ] 
that the fluid is incompressible on the t2 timescale, 



dt2P = 0. 



(72) 



so the fluid has reached the incompressible limit on the 
ti time scale. The momentum equation can be written 
as 



j„ + d} |n™/ + 1 (n*^ - n„^) | = o, 



(73) 



where we can use Eq. to substitute the velocity gra- 
dients for n*^ — Tiap- This is the usual lattice correction 
to the viscous momentum flux The kinetic modes 
make no contribution to the hydrodynamic variables, p 
and j, at long times. 

The non-equilibrium stress can be calculated by com- 
bining the stress update rule, Eqs. ((54)l and (|55p. with 
Eq. ([55]) . For example, from Eq. 



(74) 
(75) 



and from Eq. 

^*xy " Hxy = hpcl [dluy + dlu^) . 
Eliminating 11*^^ from these two equations, 

(1 - rs)K7 + hpcl {dluy + = i?.,. 



(76) 



In the general case, we again decompose the stress into 
its trace and traceless parts. 



hpcl 



1-7. 
2 



(^dlup + dluo. - '^dl^u^Sc,^ (77) 



dlfU'Sap I - Qa0 



where the random stress tensor on the macroscopic level 
is 



1 



U0 



-R. 



1 



a 13 



1 - 76 d 



(78) 



Equation (|73p can now be rewritten in terms of the vis- 
cous and fluctuating stresses 



dt2ja = dlQ 



(79) 



+5^ 



The deterministic part of the stress tensor has the desired 
Newtonian form 5] , with the usual expressions Q for the 
shear viscosity r/ and bulk viscosity C: 



c 



hpcl 1 + 7s 

2 1 
hpcl 1 



7s 
■76 



d I 



lb 



(80) 
(81) 



Combining the momentum transport on the ti and t2 
time scales we obtain the equations of fluctuating hydro- 
dynamics 0, 



dtp + da (pua) = 0, 



dt (pUa) + dp (pUaUp) + cldaP = dpQa/S 



(82) 



(83) 



C - ) djUjSap 



with random stresses Qap ■ These are Gaussian variables 
with zero mean and a covariance matrix that can be 
calculated from the analogous result on the microscopic 
level, Eq. 



{QapQ-yS) 

2mpcl 



(84) 



b'^h 



c 



2^ 

d 



This is the discrete analogue of the covariance matrix of 
the fluctuating stresses given by Landau and Lifshitz 0] . 
The delta functions in space and time that appear in 
the continuum theory are here converted into factors fe"** 
and h^^. Thus the stress fluctuations depend on the dis- 
cretization of space and time. Equation (|84p can be made 
consistent with the amplitude of fluctuating stresses in 
Ref. [H, by choosing 



m„c. 



(85) 



This is exactly the relation expected from the equation 
of state of an isothermal, ideal gas. In other words, our 
results are simultaneously consistent with macroscopic 
thermodynamics and fluctuating hydrodynamics. 



VI. CHOICE OF PARAMETERS 

The fluctuating LB model has been used to simulate 
a ran ge o f soft-matter physics, such as colloidal suspen- 
sions f33'| and polymer solutions [s^, HH] . In such cases 
it is necessary to match the LB parameters to the mass 
density, temperature, and viscosity of the molecular sys- 
tem. In addition there are two parameters that control 
the accuracy of the LB simulation without affecting the 
physics being simulated; namely the grid spacing, 6, and 
the time step, h. The grid spacing must be related to the 
characteristic length scale of the physical system. For 
example, in coupling the LB fluid to soft matter, like 
polymer chains, colloidal particles, or membranes, the 
length would be the size of the object. For flow in com- 
plex geometries, it would be the channel width, while 
for the simulations of turbulent flow, it would be the 
Kolmogorov length. This length scale, plus the desired 
spatial resolution, fixes the lattice spacing b in absolute 
units. Choosing a suitable time step then automatically 
sets the speed of sound = Csb/h, where Cs is a di- 
mensionless property of the LB model; for example, in 
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the D2Q9 and D3Q19 models = ^/iji. Typically, 
the sound speed will be unrealistically small for a dense 
liquid; however, this is not crucial since the LB method 
only runs in flow regimes where density fluctuations arc 
negligible. 

Once the length and time scales have been set, we can 
match the shear and bulk viscosities to the molecular 
system. Eqs. (|80l) and ([5T|) suggest using h and h to 
compute nondimensional viscosities from the reference 
values. 



77 



c = 



■q 



hpcl 

c 



rjh 
Qh 



hpcl pb'^t 



(86) 
(87) 



The parameters 7^ and 7^ arc then set by r) and C: 
_ 2?7- 1 

dC-1 

dC+l 



lb = 



(89) 



Small time steps therefore imply that the LB simula- 
tion is run in the over-relaxation regime. The relaxation 
rates of the kinetic modes can be chosen for convenience 
ilk = 0) or to improve the accuracy of the boundary 
conditions [2^, [2^ . 

The remaining LB parameter is the particle mass, rrip, 
which must be fixed, for a given b and h, so that the fluc- 
tuations in the LB fluid are consistent with the temper- 
ature, Eq. (|55|) . The parameter p. = rup/b^ determines 
the variance in the fluctuations [Eq. ((32|)]. 



M = 



(90) 



from which we see that too fine a grid or too large a 
time step will cause an unacceptably high noise level. A 
stable simulation will require that the time step scales as 
h cx b'^^'^'^^ or 6^/^ in three dimensions, which is slightly 
more stringent than the usual diffusive scaling, h oc b^. 



VII. CONCLUSIONS 

For models of the D3Qf 9 type, our analysis has shown 
that a fluctuating LB equation can be developed from 
statistical mechanical considerations. We have shown 
that the fluctuations are governed by the degree of 
coarse-graining, and that the relevant parameter is the 
mass of the LB particle, nip, which, for a given tempera- 
ture, is determined by the discretization of space, b, and 
time, h. The temperature appearing in the equation of 
state is identical to that which controls the fluctuations, 
as it should be. 

The beauty of the present approach is that one only 
needs to take care that the statistical properties are cor- 
rect at the LB level. The correct fluctuation-dissipation 



theorem at the Navier-Stokes level is then an automatic 
consequence of the microscopic physics. We have intro- 
duced the principle of detailed balance into the LB model, 
which is the microscopic counterpart of the fluctuation- 
dissipation theorem used in previous work 0, fiot . 
We have demonstrated all non-conserved modes must be 
thermalized in order to satisfy detailed balance' ear- 
lier implementations of the fluctuating LB model do 
not satisfy detailed balance. On the other hand, all these 
methods have been shown to be correct in the hydrody- 
namic limit. Only the stress fluctuations survive to long 
times, and the kinetic mode fluctuations become asymp- 
totically irrelevant. Nevertheless, practical simulations 
rarely probe the asymptotic limit, and then a procedure 
which is statistically correct on all length scales is clearly 
preferable. 



APPENDIX A: CONSTRAINED 
DISTRIBUTIONS 

Let us consider a constrained probability distribution 
P {{i^i}) of the following general form: 

PiW,}) (X cxp(5({z.J)) (Af) 

where 5 is a function of {t'i}, and and qj are con- 
stants. The constraints can be eliminated by making 
use of the Fourier representation of the delta function, 
S (x) = (27r) ^ J exp (ikx) dk: 

where 

S {{vi} , {kj}) = S + i^kj ViU.j - . 

(A3) 

Now, let v'f\ fej"-* denote the saddle point of S", which 
can be found by solving the linear system of equations: 



dS dS 

= ^ h i 2^ kjUi-j = 0, (A4) 



dS 



di 



(A5) 



The solution vl^"^ satisfies the constraints in Eq. (|Af p , and 
is identical to the one obtained by minimizing S, taking 
into account the constraints via Lagrange multipliers, ikj. 
The second-order Taylor expansion of S around the 
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saddle point is 



Re-introducing delta functions, we obtain the final result 



il 

where we have introduced the abbreviations 



I3ii5vi5vi + i ^ aijSi^iSkj, P {{I'l}) cx exp I ^ PaSviSn J (5 I ^ aijSvi j 

il ij \ il / j \ i / 



ft, ^ I 



2 dvidvi 



5v, 



Vi - VI 



(0) 



- — kj kj 



(A7) 

(A8) 
(A9) 



The probability distribution for Svi is then approximated 
by a Gaussian. 

The expansion of S is now inserted into Eq. (|A2|1 . Ig- 
noring the constant term, which can be absorbed in the 
normalization of P, and transforming to the new vari- 
ables 6k j, we find 

P({j/i}) cx J d{5kj)exp ^Skj'Y^aijSiyi^ 

exp PuSviSvi^ . 



(AlO) 



(All) 

Assuming the coefficients (3ij form a negative-definite 
matrix (otherwise the Gaussian approximation would not 
make sense), the saddle point is a maximum in P. 
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